function I = CalcCurrent(V, Rn, Delta, T)
    if iscolumn(V)
        V = V';
        flipOutput = 1;
    else
        flipOutput = 0;
    end
        
    e = 1.6021766e-19;
    
    %numerical integration limits
    a = 0.5;
    %a = 0.5;
    ERes = 0.00025;
    
    %approximate the DoS at important points on this mesh.
    if Delta > 0.0005
        E = [0:ERes*50:Delta/2, (Delta/2+ERes):ERes:(Delta-ERes), (Delta-0.9*ERes):(0.1*ERes):(Delta+0.9*ERes), (Delta + ERes):ERes:(Delta*2), (Delta*2+ERes) : ERes*100 : a, 5] * e;
    else
        res = Delta / 500
        E = [0:res:3*Delta, 3*Delta+ERes : ERes * 100 : a, 5] * e;
    end
    E(E == Delta * e) = [];
    
    %  E = [0:1*ERes:1*e];
    dosCalc(1, :) = E;
    dosCalc(2, :) = DoS(E, Delta * e);
    
    
    %now interpolate the mesh for actual points. 
    factor = 1;
    E = (-a : ERes/factor: a) * e;
    %from the calculation, we interpolate the values we will use to do the calculation.   
    ETable = ones(length(E), length(V));
    for i = 1:length(V)
        ETable(:, i) = E + e*V(i);
    end
    
    dosTable = interp1(dosCalc(1, :), dosCalc(2, :), abs(ETable), 'linear');
    %dosTable = DoS(ETable, Delta);
    occTable = GetOccupation(ETable, T);
    dosTableV0 = repmat(interp1(dosCalc(1, :), dosCalc(2, :), abs(E))', 1, length(V));
    occTableV0 = repmat(GetOccupation(E, T)', 1, length(V));
    I = trapz((dosTableV0 .* dosTable .* (occTableV0 - occTable)) * ERes*e/factor, 1) / (e * Rn);
    
    if flipOutput
        I = I';
    end
end

function I = CalcCurrentAuto(V, Rn, Delta, T)
    if iscolumn(V)
        V = V';
        flipOutput = 1;
    else
        flipOutput = 0;
    end
        
    e = 1.6021766e-19;
    
    Delta = Delta * e;   
    
    I = integral(@(E) DoS(E, Delta) .* DoS(E+e*V, Delta) .* (GetOccupation(E, T) - GetOccupation(E+e*V, T)), -Inf, Inf) / (e*Rn);
end

function N = DoS(E, Delta)
    %assume E is monotonically increasing.
    idx = find(abs(E) > Delta);
    otherIdx = find(abs(E) <= Delta);
        
    N = zeros(size(E));
   % N = mask * 2*E(idxA:idxB) .* ellipticK(E(idxA:idxB) / Delta) / (pi*Delta);;
    
    N(idx) = 2 * ellipticK((Delta./E(idx)).^2) / pi;
    N(otherIdx) = 2*abs(E(otherIdx)) .* ellipticK((E(otherIdx) / Delta).^2) / (pi*Delta);
       
     N((N == Inf)) = 5;
end

% function N = DoS(E, Delta)
%     N = ones(1, length(E));
% end

